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Flows of hard granular materials depend strongly on the interparticle friction coefficient /ip and 
on the inertial number T, which characterizes proximity to the jamming transition where flow stops. 

Guided by numerical simulations, we derive the phase diagram of dense inertial flow of spherical 
particles, finding three regimes for 10“^ 2i < 10“^: frictionless, frictional sliding, and rolling. 

These are distinguished by the dominant means of energy dissipation, changing from collisional to 
sliding friction, and back to collisional, as /tp increases from zero at constant T. The three regimes 
differ in their kinetics and rheology; in particular, the velocity fluctuations and the stress ratio both 
display non-monotonic behavior with /tp, corresponding to transitions between the three regimes 
of flow. We rationalize the phase boundaries between these regimes, show that energy balance 
yields scaling relations between microscopic properties in each of them, and derive the strain scale 
at which particles lose memory of their velocity. For the frictional sliding regime most relevant 
experimentally, we find for I > 10“^'® that the growth of the macroscopic friction with I is 
induced by an increase of collisional dissipation. This implies in that range that /i(X) —/i(0) ~ 
where b ~ 0.2 is an exponent that characterizes both the dimensionless velocity fluctuations C ~ 
and the density of sliding contacts y ~ T**. 


Dense flows of granular media are central to many in¬ 
dustrial processes and geophysical phenomena, includ¬ 
ing landslides and earthquakes m- At a fundamental 
level, describing such driven materials remains a chal¬ 
lenge, in particular near the jamming transition where 
crowding effects become dominant and flow stops. In the 
last decade, progress was made by considering the limit of 
perfectly rigid grains, for which dimensional analysis im¬ 
plies that the strain rate e, the pressure P and the grain 
density p can only affect flows via the inertial number 
I = €D^/~pJP, where D is grain diameter Hi]. In partic¬ 
ular, for stationary flows the packing fraction cj) and stress 
anisotropy /t = cr/P, where a is the shear stress, are func¬ 
tions of I. From the constitutive relations 4>(P) and pfl) 
the flow profile can be explained in simple geometries 
[1[7HS|. Here we focus on dense flows I < 0.1 for which 
the networks of contacts between grains span the system 
and particle motion is strongly correlated mm, and 
do not consider the quasi-static regime I < 10“"*’ where 
flow appears intermittent [illlMlj- In this intermediate 
range one finds 

pil) =Tc + af, , (fil) = (l)c-a^ , ( 1 ) 

where and (fc are non-universal and depend on de¬ 
tails of the grains. Experiments on glass beads and sand 
find exponents « 1, consistent with numerical 

simulations using frictional particles reporting = 0.81 
and = 0.87 [11]. Despite their importance, constitu¬ 
tive laws Eq]!] remain empirical. Building a microscopic 
framework to explain them would shed light on a range 
of debated issues, including transient phenomena mm, 
non-local effects [IMS], and the presence of S-shaped 
flow curves when particles are soft [^01422) . 

To make progress, it is natural to consider the limit¬ 
ing case where particles are frictionless, a situation that 
has received considerable attention in the jamming liter- 



FIG. 1. (Golor online) Phase diagram of dense homogeneous 
inertial frictional flow. In the frictionless and rolling regimes, 
most energy is dissipated by inelastic collisions, while in the 
frictional sliding regime energy dissipation is dominated by 
sliding. Along the phase boundary, grains dissipate equal 
amounts of energy in collisions and in sliding. For X > 0.1, 
one enters the dilute regime [9]. The dashed line has slope 2. 


ature . For hard particles, two geometrical results 

key for inertial flows are as follows. First, as the density 
increases, the network of contacts becomes more coordi¬ 
nated, implying that motion becomes more constrained. 
This leads to a divergence of the velocity fluctuations 
{6V) when constraints are sufficient to jam the material 
PTIlddj . Thus the contact network acts as a lever, whose 
amplitude is characterized by the dimensionless number 
C = {5V)/{eD). At the same time, the rate at which 
new contacts are made increases, and the creation of each 
contact affects motion on a growing length scale. These 
effects imply that velocity fluctuations decorrelate on a 
strain scale Cy that vanishes at jamming |3T]. The the¬ 
ory of Ref. |3T] , which uses the fact that dissipation can 
only occur in collisions for frictionless particles, predicts 
= 0.35, C ~ and ty ^ X. Encourag¬ 

ingly, these results agree with the numerics of Ref. [32] , 
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FIG. 2. (Color online) Ratio of dissipation due to sliding, 
'Dsiid, to dissipation from collisions, 'Dcoii, vs fip. The triangle 
has slope 1. 

which found « 0.38 and C ~ How¬ 

ever, Ufi and differ significantly from their values for 
frictional grains stated above, suggesting the presence 
of different universality classes. Currently, why friction 
qualitatively affects flows and potentially leads to sev¬ 
eral universality classes, how many universality classes 
exist, and what differs between them microscopically are 
unresolved questions. 

In this work we use numerical simulations to answer 
these questions. We systematically study dense flows 
over a large range of X and fXp. By focusing on the mi¬ 
croscopic cause of dissipation, we show the existence of 
three universality classes, as illustrated in Fig. At low 
friction, there exists a frictionless regime in quantitative 
agreement with the theory of Ref. m, in particular 
we establish that ~ I. As the friction increases, one 
enters the frictional sliding regime, where dissipation is 
dominated by sliding at contacts instead of collisions, and 
for which ~ I holds true but C ~ X~^ with b « 0.22. 
We relate the exponent b to the density of sliding con¬ 
tacts, X ^ 21^. Most importantly, we show that although 
the value of in Eq{^ reflects sliding dissipation, the 
dependence of fi with X is governed by collisional dissi¬ 
pation when X > 10“^ ®, leading to = 1 — 2b. Fi¬ 
nally, at even larger fj.p one enters a rolling regime where 
dissipation is once again dominated by collisions, and 
where exponents are consistent with those of frictionless 
particles, both for kinetic observables and constitutive 
laws. We derive the phase boundary between the fric¬ 
tionless and the frictional sliding regime. Overall, our 
work explains why friction qualitatively changes physical 
properties, and paves the way for a future comprehensive 
microscopic theory of dense granular flows. 

Numerical Protocol — To model inertial flow of fric¬ 
tional particles, we use a standard discrete element 
method |33j in two dimensions, described in more de¬ 
tail in Appendix 1. Collisions are computed by modeling 
grains as stiff viscoelastic disks: when grains overlap at a 
contact a, they experience elastic and viscous forces /®, 
and /^, respectively, leading to a restitution coefficient 
which we choose to be e = 0.1 [53]. The tangential (nor¬ 



FIG. 3. (Color online) Lever C vs I, (a) for < 0.3, and (b) 
for /ip > 0.3. The dashed lines are oc , while the dotted 

line is oc 

mal) components fd’" (fd) are restricted by Coulomb 
friction to satisfy Ifd^l < t^pfdt contacts that saturate 
this constraint are said to be sliding, while those that 
obey a strict inequality are said to be rolling. 

Shear is imposed with rough walls bounding the up¬ 
per and lower edges on an x—periodic domain. We per¬ 
form our numerics at imposed global shear rate and con¬ 
stant pressure, following a system preparation described 
in Appendix 1. We discard data that do not satisfy 
strict criteria for homogeneity of the flow, as specified in 
Appendix 1. Grain stiffness is such that relative defor¬ 
mation at contacts is A « 10“® ®, within the rigid limit 
established previously [5|, and system size is large enough 
to ensure the absence of finite-size effects. Independence 
of our results with respect to A, e, and N is shown in 
Appendix 2. 

Partitioning dissipated power — Frictional particles can 
dissipate energy either through inelastic collisions, at a 
rate Pcoih or by sliding at frictional contacts, at a rate 
Vsiid. In our contact model, inelasticity is due to the 
viscous component of contact forces; therefore the colli¬ 
sional dissipation rate I?coH can be written 

^ Y. + E ’ (2) 

a^C aGCfi 

where Ua is the relative velocity at contact a, decom¬ 
posed into normal and tangential components, C/^ and 
HJ. Here C denotes all contacts, of number and 
Cr denotes rolling contacts. The dissipation rate due to 
sliding is 

Vsud =Y ’ (3) 

ctGCs 

where Cs is the set of sliding contacts. In steady state, 
dissipation must balance the work done at the boundaries 
[35] . The energy input from the shear stress is VLae, where 
H is the system volume, and for large systems, additional 
contributions from fluctuations of the normal position 
of the wall are insignificant. We define dimensionless 
dissipation rates per particle VcoU = Vcou/i^pi), X>sUd = 
Vsiid/i^pe), so that [55] 

P = 'Xt coll pit slid- ( 4 ) 
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FIG. 4. (Color online) Autocorrelation of particle velocities, 
C{t) = {V^{Q)v^{e))/{V^{Qf) for /Tp = 0.3 and indicated 
T. (b) C{t) vs e/e„. In (b), unfilled symbols correspond to 
strains larger than 0.01, not used for fitting, and the solid line 
shows the fitted form. 


FIG. 5. (Color online) (a) Decorrelation strain scale, vs 
I for selected jj,p. The dotted, dashed, and dot-dashed lines 
are oc and , respectively, (b) The fraction of 

sliding contacts, Xj ''S for various /ip (symbols as in Fig. 
3). Dotted, dashed, and dot-dashed lines have slopes 0.43, 
0.41, and 0.27, respectively. 


To investigate which source of dissipation dominates in 
Eq.m we consider the ratio 'DgHj^/'Dcoih shown in Figj^ 
As expected, collisional dissipation dominates in the fric¬ 
tionless limit, but sliding dissipation becomes more im¬ 
portant as /ip is increased, and becomes dominant at in¬ 
termediate friction coefficients and small inertial num¬ 
ber, consistent with earlier simulations for /ip = 0.3 |36j . 
Strikingly, the dependence on /ip is non-monotonic: when 
/ip reaches « 0.2, this trend abruptly reverses, and 
'Dsiid/'Dcoii decreases with /ip, implying that collisional 
dissipation dominates as /ip —>■ oo. 

To define phase boundaries, we use the inertial number 
at which = !> resulting in the phase diagram 

of Fig.[^ From the non-monotonicity of 'DsUd/'Dcoii with 
/ip, this leads to two phase boundaries merging at I « 
0.1, where the dense flow regime ends 111136]. This defines 
three flow regimes: frictionless, frictional sliding, and 
rolling, where sliding dissipation dominates only in the 
intermediary regime. Later in this work, we will show 
that this phase diagram correctly classifies kinetics as 
well as constitutive laws. 

Connecting dissipation to key kinetic observables — In 
the rigid limit, collisions become very short in duration, 
and the power dissipated in collisions can be expressed in 
terms of microscopic observables m, as we now recall. 
Each time a particle changes its direction with respect 
to its neighbors, a finite fraction of its kinetic energy 
^ mSV'^ must be dissipated, where m is the particle mass 
(we consider finite restitution e < 1). Since Cy is the 
characteristic strain at which velocities decorrelate, this 
occurs at a rate cx e/e„, thus Vcoii ^ N{e/ey)m{5V'^) and 


T^coii oc 


N{e/ey)mSV'^ 

ND‘^pe 


oc - 




(5) 


The rate of sliding dissipation can be directly esti¬ 
mated from its microscopic expression, Eq. (|^. We 
assume that the force at the sliding contact is typical, 
i.e. I/Jl = fJ^pfa ~ and that the sliding 

velocity is of the order of the velocity fluctuation, i.e. 
|C/J| ^ 6V. These assumptions hold true in the slid¬ 
ing frictional regime where they matter (they eventually 


break down in the rolling regime where sliding contacts 
become rare and atypical, see Appendix 3). We get the 
estimate 


NM\f^\)s{\U^\)s 

ND'^pe 


'^slid OC 




( 6 ) 


where {■)s denotes an average over sliding contacts, 
whose fraction is y. Using Eqs. ( 4|5|6 ) we now get the 
following constraints on the different regimes: 


/t ~ ^ Frictionless, Rolling (7) 

/t ~ pi-pX^^ Frictional Sliding (8) 


We now test these scaling relations and use them to com¬ 
pute the boundary of the frictionless regime. 

Measuring kinetic observables- We measure the lever 
effect defined as £ = {6V)/{eD), where {6V) is the typ¬ 
ical magnitude of velocity fluctuation about the mean 
velocity profile [37]. Our results are shown in Figj^ For 
any /ip, C grows as I —>■ 0. In the frictionless limit, we 
find C (X in agreement with earlier results [35] 

and the prediction m- A striking result is that the am¬ 
plitude of this growth is non-monotonic in /tp, with a 
minimum around /ip « 0.2, thus closely paralleling the 
phase diagram of Figj^ Moreover, in the /ip —)■ oo limit, 
the divergence is again close to £ cx In contrast, 

curves that are fully in the frictional sliding regime, as oc¬ 
curs for /ip = 0.1 or /ip = 0.3, are well fitted by £ cx 
with b = 0.22, close to experiments finding £ cx X~^/^ 

[niiH]. 

We now turn to the strain scale Cy beyond which a 
particle loses memory of its velocity. It can be ex¬ 
tracted from the decay of the autocorrelation function 
[33] C{e) = (I4^(0)I^^(e)), where we use the vertical com¬ 
ponent of velocity at particle i, U/, averaged over all 
particles and initial time steps. The normalized correla¬ 
tion function C{e) = C'(e)/C'(0) is shown for pLp = 0.02 
and various X in Fig. We see that beyond a scale Cy, 
C(e) decays as a power-law, as observed numerically in 
over-damped suspensions [SHIM]. For all /ip, C(e) has 
a similar form; we find that for e < 10“^ it is well-fitted 
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by [1 + {e/ey{I)Y)\ with v = 1.1 and rj dependent 
on fip. By rescaling e to obtain a collapse, shown in 
Fig. we obtain the scale e„. Repeating this process 
for all Up leads to the results shown in Fig. for selected 
fip', other Up are shown in Appendix 4. observe that 
for all T and all fip, we have approximately ey « I, al¬ 
though our best exponent for the rolling regime is closer 
to €v « This result is thus in excellent agreement 

with the prediction of [3T] for the frictionless case, and 
also with experimental measurements finding e^, ~ X |38j . 
We recall below our previous argument, which we expect 
to hold more generally for frictional particles. 

Finally, the fraction of sliding contacts x is shown in 
Fig.j^D. For each fip, x decays as X —>• 0. In the friction¬ 
less regime, y « 1 , as expected, while in the frictional 
sliding and rolling regimes, y decays as a power-law as X 
is decreased. For the frictional regime, such as fip = 0.3, 
data are well-fitted by y ^ while in the rolling 

regime we find a sharper decay, y ^ X^'^® for fip = 10 . 

Our results on microscopic quantities are summarized 
in Table |lj We see that the scaling relations Eqs. 
are consistent with the data. 

Regime boundaries - We can now estimate when the 
frictionless regime breaks down. Since in that regime 
C ~ X“° ®, e„ ~ X, and y ^ 1, we have according to 


Eqs. (Sjei Vsiid/'Dcoii ^ ftpX consistent with Eigj2j 
The frictionless regime must break down at an inertial 
number X^ where this ratio is of order one, yielding X^ ^ 
fip in agreement with Eig{^ 

Inside the sliding regime, we have « X. To de¬ 
termine the transition to a rolling regime, we note from 
Eq. ® that Vsiid/Vcoii II{ILY We ob¬ 

serve chat the product y/ip decays with large fip at fixed 
X (data not shown). Thus, although the dissipation of 
each sliding contact grows with ftp, fewer and fewer con¬ 
tacts slide as fip becomes very large, and the latter effect 
dominates when fip is large enough. This qualitatively 
explains the observed non-monotonic behavior with fip. 

Constitutive Relations- Experimentally, the most ac- 


Regime 

Relation 

Prediction 

Measured 

Frictionless 

£ ~X-'’ 

b= 1/2 

c = 1 

b = 0.50 

c= 1.10 

Frictional 

Sliding 

£ ~X-'’ 

y~x" 

Sfi ~ 

c = 1 

d = b 

ctp = l — 2b 

b = 0.22 

c = 0.95 

d = 0.27 

Up = 0.6 

Rolling 

y~x" 

c = 1 

b = 0.50 

c= 1.28 

d = 0.43 


TABLE I. Summary of scaling behavior. Predictions for the 
frictionless regime are quoted from the theory of |31| . while 
the other predictions are Eqs. (TjSl. In the frictional sliding 
regime, scalings are taken for the extremal value fip = 0.3, 
while for the rolling regime, scalings are taken from fip = 10. 


cessible quantities are the constitutive relations /i(X) and 
4>{I), which we show in Eig. To discuss universality 
classes, it would seem appropriate to measure the ex¬ 
ponents Op and Q !0 entering Eq. Q. However, these 
exponents are much harder to measure than those sum¬ 
marized in Table I, because of finite-size effects in the 
fitting parameters fiy and (jic |32j . Instead, we simply 
consider the cases fip = 0, fip = 0.3, and fip = 10 for 
which our data are respectively in the frictionless, fric¬ 
tional sliding, and rolling regimes. In the inset to Figj^ 
we show that Sfiil) = /i(X) — fi^ is nearly identical in the 
frictionless and rolling regimes (the points overlap), and 
definitely distinct from its behavior in the frictional slid¬ 
ing regime. This observation supports further our claim 
for three distinct universality classes. 

We now argue that in Erictional Sliding regime, the 
exponent describing the evolution of the macroscopic 
friction with inertial number as defined in Eq.Q can be 
deduced from the exponent b characterizing the velocity 
fluctuations. From Eq. the partition of dissipation 
is also a partition of fi. As shown in Eig|^, we observe 
that the contribution from sliding is nearly independent 
of X, while the contribution from collisions is vanishing as 
X —>■ 0. So long as the variation in sliding dissipation with 
X is negligible, this implies that p.(X) — fi^ is dominated 
by collisional dissipation, even in the frictional sliding 
regime. We find that this is the case for X > 10“^ ® 
(data not shown). This facilitates precise measurement 
of ^(X) — fic in this range, with which we obtain the 
meas urem ent = 0.6 for fip = 0.3. Moreover, using 
Eqs. (^) we obtain 


Sfi{I) oc-~ XX 


2 _ jl-'ib 


Frictional Sliding, 


(9) 


implying = 1 — 2b. Our prediction « 0.6 is in 
reasonable agreement with the previous measurement of 
m where « 0 . 8 , considering the restricted range of 
inertial number where we expect this power-law behavior 
to hold. 

Scaling argument for the characteristic strain scale- 
The relation « X can be rationalized by a generaliza¬ 
tion of the argument in m- We use the geometrical fact 
that in dense flows, when a grain has an unbalanced net 
force, F, the ensuing motion will tend to make the re¬ 
maining contacts of the grain align along F (see Fig0. 
Since forces are repulsive, this further increases the un¬ 
balanced force and further accelerates the grain. The in¬ 
crease in net force is proportional to the typical contact 
force, , as well as to the rotation of the contacts, 

of magnitude ~ Cde, thus 

( 10 ) 

where C is the dimensionless magnitude of the velocity 
fluctuation. This equation can also be derived formally, 
as shown in Appendix 5. In inertial flow, unbalanced 
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= 0'02 


▼ '^slidy Mp ~ 0.3 


l^-p = 0.02 

(c) 

-m-Vcall, lip = 0.3 


FIG. 6. (Color online) (a) Volume fraction 0 vs I for various fip (symbols as in Fig. 3). (b) Effective friction fj, vs I (symbols 
as in Fig. 3). Inset shows non-monotonic behavior of fi{T) — jj,c, for = 0,0.3,10. (c) Decomposition of fj, into collisional and 
frictional components. Both lines have slope 0.6. 


forces are proportional to accelerations, F = pl'^dC/de, 
which leads to 


d^C 

de^ 


oc 


£ 

12 - 


( 11 ) 


Eq[n] indicates that there is a characteristic strain scale 
tp « I in which a velocity fluctuation grows by an 
amount proportional to its initial magnitude. In steady 
flow, such growth must be destroyed by collisions on the 
same strain scale, since the latter reorganize the direc¬ 
tion of particle motion. Hence this is indeed the scale 
of decorrelation of particle velocities. (At very large p,p, 
the direction of a contact force can be unrelated to the 
contact direction, and corrections to our argument are 
plausible.) 

Discussion- In this work we have shown that dense in¬ 
ertial granular flows can be classified into three regimes, 
in a phase diagram spanned by the friction coefficient /ip 
and the distance to jamming, characterized by the iner¬ 
tial number I. By considering the microscopic cause of 
dissipation, we have shown that its nature must change 
as the friction coefhcient /tp increases from zero. One 
eventually leaves the frictionless regime to enter in the 



FIG. 7. (Golor online) Illustration of geometrical nonlinearity. 
If the central grain has an unbalanced force as indicated by 
the arrow, then the ensuing flow will tend to align the contact 
normals of the dominant contact forces (thick lines), i.e. the 
angle (f) will increase. Geometrically, d(f>/dt oc V, the velocity 
of the particle. 


frictional sliding regime, where both the kinetics and con¬ 
stitutive relations differ. As /ip increases further, fewer 
contacts slip, and one enters the rolling regime where col¬ 
lisions once again dominate dissipation, and where expo¬ 
nents are consistent with that of the frictionless regime. 

Experimentally, these results could be tested by mea¬ 
suring the correlation function C{e) = {V^{0)V^(e)), 
which captures both the lever amplitude C (at e = 0) 
and the strain scale e„. This will require a sufficient res¬ 
olution in the strain e that can be probed. Varying the 
friction coefficient in these studies would also be valuable. 

On the theoretical level, a complete theory of the fric¬ 
tional sliding regime, the most important in practice, is 
still lacking. Here we have proposed a scaling descrip¬ 
tion relating the singularities in the constitutive law /i(I) 
to those in the kinetic observables CpiT), C{X) and x(X), 
which can all be expressed in terms of a single unknown 
exponent b. A key challenge for the future is to predict 
the value of b. Moreover, our arguments are mean-field in 
nature, as they assume that dissipation occurs rather ho¬ 
mogeneously in space, and that velocity fluctuations are 
described by a single scale C. Although there is evidence 
that such mean-field arguments are exact for frictionless 
particles |31| , they may be only approximate in the fric¬ 
tional case where intermittent strain localization is some¬ 
times reported 0 M- Concerning the rolling regime, 
why it has the same scaling exponents as the frictionless 
regime also needs to be clarified further, beyond their 
similarity in dissipation mechanism established here. 

Finally, this work could be extended in several direc¬ 
tions. It would be very interesting to measure the kinetic 
quantities presented here in the intermittent quasi-static 
regime of very slow flows I < 10’ [3 [MU]- Similar 

extensions could be done with respect to particle shape, 
where local ordering is important HOI m], and parti¬ 
cle softness, where the flow curve can become sigmoidal, 
leading to hysteresis [201122] . Last, over-damped suspen¬ 
sions present the same problem as inertial flows: various 
numerical studies have focused on frictionless particles 
[muziiiiHis], which appear consistent with the theory 
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developed in m- Together with Eq. the theory pre¬ 
dicts that frictional sliding should dominate over viscous 
dissipation when r]Qi/P /i^, where 770 is the viscosity 
of the solvent. It is currently unclear whether this tran¬ 
sition qualitatively affects physical properties, as experi¬ 
ments [46l |47] and numerics [48] with friction are reason¬ 
ably compatible with the frictionless theory. Numerically 
building a phase diagram analogous to Fig. compar¬ 
ing the amplitude of sliding dissipation to other sources, 
would resolve this issue. 
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I. APPENDIX 

In these Appendices, we provide additional details of 
our results. Appendix 1 describes details of the numer¬ 
ical simulations, while Appendix 2 shows that our main 
results are independent of the grain stiffness, restitution 
coefficient, system size, and numerical integration time- 
step. Appendix 3 shows atypical behavior of sliding ve¬ 
locity and sliding force in the rolling regime. Appendix 
4 shows the velocity autocorrelation function for several 
values of fj,p and e„ for all values of /ip considered. Ap¬ 
pendix 5 computes the effect of geometrical nonlinearity 
during flow. 


Appendix 1. Numerical Simulations 

Simulations are performed with a standard Discrete El¬ 
ement Method code [35], which integrates Newton’s equa¬ 
tions of motion for each grain with Verlet time-stepping. 
We focus on two dimensions, as empirically exponents do 
not appear to depend on dimension; see |31| for a review 
of the literature on this point. Collisions are computed by 
modeling grains as viscoelastic disks: when grains overlap 
at a contact a, they experience elastic and viscous forces 
/®, and /^. The coefficient of the viscous force is chosen 
to obtain a restitution coefficient e = 0.1 in binary colli¬ 
sions; away from the singular limit e —>■ 1 that we do not 
consider, varying this coefficient does not strongly affect 
our results, as shown in Section 2. These forces can be 
decomposed into their contributions normal to the con¬ 
tact, and and tangential to the contact, 
and ja''’. The tangential force is imposed to stay inside 


the Coulomb cone, |/a^| < fJ-pfa ■ Contacts that saturate 
the Coulomb constraint are said to be sliding. 

The grains are polydisperse with equal numbers of di¬ 
ameter (0.82,0.94,1.06,1.18), the same mixture used in 
[33] . Previous work established that the polydispersity 
does not affect /ic in simple shear flow, even over a huge 
range of polydispersity |50j . Results in the main text 
are reported for a value of grain stiffness such that the 
grain relative deformation is set to A = 10“^ ®, ap¬ 
propriate for some materials [3]. This is well within 
the range A < 10“® in which rheological results are 
independent of A, as previously established in simula¬ 
tions of inertial flow of frictionless and frictional particles 
[SKUIIITIIMIIIH]. This result is verified by the explicit de¬ 
pendence of our phase diagram on A, reported in Section 
2 for A e {10-®-®, 10-2-8,10-1-8}. We studied three sys¬ 
tem sizes N e {1000,1800,3700}. Results are reported 
for the largest N; the absence of finite-size effects is es¬ 
tablished in Section 2. The shortest time-sca le of the 
dynamics is the microscopic elastic time-scale y/mjk; we 
chose our numerical time-step At such that it never ex¬ 
ceeds Q.Q&^/rnJk. This ensures that binary collisions are 
resolved with > 15 steps, and the much slower multi¬ 
body collisions typical of dense flow will be resolved in 
even greater detail. Independence of our results with re¬ 
spect to At is shown in Appendix 2. 

The square domain of size x Ly is periodic in the 
x-direction and has upper and lower walls. The walls 
are created from the same polydisperse mixture as the 
bulk, staggered to create roughness. The walls obey an 
equation of motion 


M 


dt'^ 


dr _ ^ 

Vi, — ^bulk 

at 


( 12 ) 


where M is mass, r/ is a damping coefficient, is the 
force from the bulk of the packing, and F^xt is an external 
applied force. The bulk-wall interactions are via contact 
forces, exactly as in the bulk. The external force in the y- 
direction is constant, such that on the top 

(j-) and bottom (-) walls. In the x-direction, the external 
force is chosen to impose a constant velocity ±t4j, and 
hence a constant global shear rate e = 214,/Lj,, up to 
fluctuations in Ly. 

We seek to make the flow as ho mogeneous as possible. 
Following [^, we set rj = •\/mk, where k is the spring 
constant for particle-particle elastic interactions, and m 
the mean particle mass. We tested the dependence of the 
results on M. When M/m ~ 1, the wall equation Eq.( 12) 
is dominated by the viscous term, and can exhibit long 
transients. We therefore set M/m = 50, so that the 
wall density and particle density are the same order; this 
minimized transients. 


With this choice of wall parameters, we find that 
steady states are achieved where the relative pressure 
fluctuations range from 1% at I ^ IQ-® to 20% at 
I ~ 0.1; thus the mean particle overlap A cx P/k is 
fixed to within this precision. 








FIG. 8. (Color online) A^-dependence of observables, (a) Dissipation ratio vs X. (b) x vs X. (c) C, vs X. 





FIG. 9. (Color online) e-dependence of observables, (a) Dissipation ratio vs X. (b) x vs X. (c) L vs X. 


X 


10-1 

10-2 

10-3 

10 -^ 



FIG. 10. (Color online) (a) A-dependence of phase diagram for A = 10 ® ® (squares, solid), 10 (diamonds, dashed), and 
10 - 1 .s (triangles, dash-dotted), (b) Absence of dependence of C, on time step At for Hp — 0.02. 


To prepare homogeneous steady states, initially 
isotropic packings are created from a gas at volume frac¬ 
tion (f>Q = 0.5, and then sheared for a pre-strain eg. As 
discussed below, an analysis of Eq{T^ leads to our choos¬ 
ing eg = max(0.2, A“^/^X), which we checked ensures 
that a steady state is reached. After this initial strain, 
without collecting data, we strain the systems for e = 0.3, 
collecting data every (5e = 5 x 10“"^. 

In all cases, we discard runs that are not sufficiently 
homogeneous. As a first criterion, we exclude simulation 
runs where the mean velocity profile has a shear band. 


As a second criterion, we find for certain parameter val¬ 
ues that resonant elastic waves bounce back and forth 
between the walls at very high frequency, as discussed in 
|51) . Resolution of these waves requires a much smaller 
time step than is needed otherwise, so we do not include 
these runs. Details of these criteria follow. 

To determine an appropriate pre-strain scale eg, con¬ 
sider the j/-direction bulk-wall force on on the top wall, 
^buik- "This is a spring-like force, since it results from 
the elastic interactions between the particles adjacent 
to the wall, and the wall itself, but with a nontriv- 
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ial spring constant. It can be estimated from the law 
(j)c — fp ^ 1-°‘* ■ Indeed, linearizing this law around a 
mean volume fraction (j) and mean pressure P, we find 
P - P (X - Lyl2)lLy, where 

Ly is the mean thickness of the domain and we used 
{(j) — = —{Vw — Lyl2)l{Lyl2). Hence the bulk-wall 

force is approximately F^Jf. ~ ~ -kwiVu, - Ly 12) 

with ku, ~ i'^{4>c — ~ P{,(t>c — $)~^■ The strain 

scale associated with the damping term in Eq. is then 
eo ~ iil/kw — $), where A = p/k. We con¬ 

servatively take eo = max(0.2, A~^/^I). 

Our two criteria for ensuring homogeneity of the flows 
are that there is no static shear band, and that the walls 
are not in resonant motion. To test for a shear band, 
we compute the deviation of the mean velocity profile 
from a linear one, 5v{y) = v{y) — ey, and compute its 
normalized standard deviation, {{6v{y) — {6v))‘^)y/{Lye)^. 
For a perfect shear band, this is l/-\/l2 = 0.29; we discard 
runs where it exceeds 0.2. 

For certain parameter values, resonant elastic waves 
bounce back and forth between the walls at very high fre¬ 
quency, as discussed in m- Resolution of these waves re¬ 
quires a much smaller time step than is needed otherwise, 
and in our code they display an unphysical alternation of 
the velocity of the wall from positive to negative values 
at each strain increment where we save data. Therefore 
we compute a normalized numerical derivative of the ver¬ 
tical wall velocity, O = (t4j(e + Ae) — I4,(e))/(Ae((5F)), 
where V^, is the wall velocity (for brevity, here we include 
only one wall), and (SV) is the velocity scale of grains in 
the bulk, computed from their fluctuations. We find that 
for well-behaved runs, OP ^ 1, while for numerically un¬ 
stable ones, OX > 10000. Therefore we exclude runs with 
OX > 2000. We checked that the few runs so excluded 
agree in their location in parameter space with the theory 
of [5T] . 


Appendix 2. Dependence of results on A,e,N, and At 

In the main text, we reported results for grain rel¬ 
ative deformation A « 10“^’®, number of particles 
N « 3700, r estitu tion coefficient e = 0.1, and time step 
At « O.OO^/m/fc. Here we discuss how our results de¬ 
pend on these choices. We show representative plots for 
N « 1000,1800,3700 at several values of pp and sev¬ 
eral quantities in Figj^ We see that x and Vsiid/X>coii 
are independent of N at the two largest values stud¬ 
ied. Since velocity fluctuations are suppressed at the 
wall, L displays an expected mild, systematic depen¬ 
dence. This behavior is representative for all values of 
Pp. In all cases the minor dependence in C does not affect 
the scaling behavior of observables, and in particular the 
phase diagram is not affected. Similarly, representative 


plots for e = 0.5 in small systems N « 1000 show that 
Pslid/Pcolli nnd L display only a weak dependence on 
the restitution coefficient. 

Although the grain relative deformation A « 10“® ® 
is well within the rigid limit established in previous 
work [3 na [HI [321 ss], we checked how our phase 
diagram depends on A G {10“®®, 10“^®, 10“® ®}, as 
shown in Fig{T§i. For the two smallest values of A, 
the frictionless-frictional-sliding transition is indepen¬ 
dent of this value above the quasistatic regime, and the 
frictional-sliding-rolling transition displays only a very 
small dependence. 

Finally, in FigllO}) we show independence of our results 
on the time-step At, which was halved for a set of sim¬ 
ulations with pp = 0.02, which includes both frictionless 
and frictional sliding regimes. 

Appendix 3. Sliding dissipation in rolling regime 

In the rolling regime, only a small subset of con¬ 
tacts are sliding, as shown in Fig.5b of the main text. 
This raises the possibility that the forces and veloci¬ 
ties at these contacts may be atypical of the system 
as a whole. We investigate this with the quantities 
C = and £t = (|[/^|)s/(De), where 

{■)s denotes an average over sliding contacts. As shown 
in FigfTl] atypical behavior of sliding contacts is indeed 
shown for large fip, and in fact we find in this regime that 
both C and £t j £ show power-law behavior. In partic¬ 
ular, for pp = 2 we find ( ~ I® ®® and £tI£ X“°-^®, 

while for pp = 10 we find ^ ^ I®-^® and £tI£ ~ 1“° ®'®. 
The quantities C and £t/£ would be needed for an accu¬ 
rate scaling estimate of sliding dissipation in the rolling 
regime. 

Appendix 4. Velocity autocorrelation function 

In the main text we introduced the autocorrelation 
function 

C(e) = (H/(0)H/(e)), (13) 

in terms of the vertical component of velocity at parti¬ 
cle z, H/, averaged over all particles and all initial time 
steps for a given strain increment e. In Fig |12[ we show 
C'(e)/C'(0) for several values of Pp. The values of ry are 
listed in Tab le [lH and the resulting values of e«(I) are 
shown in Fig|13| 


Pp 

0 

0.002 

0.02 

0.1 

0.3 

0.5 

0.8 

2 

10 

V 

1.2 

1.2 

0.75 

0.5 

0.5 

0.55 

0.65 

0.75 

0.75 


TABLE II. Exponent rj in fit of C(e)/C(0). 
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FIG. 11. (Color online) (a) ^ vs X. Dotted and dashed lines have slopes 0.43 and 0.31, respectively, (b) Ct C vs I. Dotted 


and dashed lines have slopes -0.34 and -0.28, respectively. 



e e 





^ip = 0.8 ^ = 10 





FIG. 12. (Color online) Autocorrelation of particle velocities, C'(e) = (Fi^(0)y/(e)) for indicated /Xp and I from 10“^'^ (down 
triangles) to 10~^'^ (circles) (symbols as in Fig 4a of the main text), (a-d) C{t)/C{Q) vs e. (e-h) C{t)/C{Q) vs e/e„. In (e-h), 
unfilled symbols correspond to strains larger than 0.01, not used for fitting, and the solid line shows the htted form. 


Appendix 5. Geometrical nonlinearity in flow 

We aim to compute how forces evolve along flow of hard particles. In particular, we consider flow along a floppy 
mode, where the relative velocity at a contact, Uij, is zero at rolling contacts, and only transverse at the sliding 
contacts. The definition of Uij, in particular the motion of j relative to i at their mutual contact, is 

'^ij — T) T ( 1) ^ 

where Awij = RjCOj + RiUJi- This holds both in d = 2 and d = 3, where the cross product in 2D is defined as 
V X u = ViCijUj in terms of the Levi-civita symbol ei 2 = —£21 = IjCii = £22 = 0. In 2D uj becomes a scalar, and 
iXxfi = u!Xn is the vector ujckini. In this way we can handle both cases simultaneously. 

By multiplying Uij along an arbitrary set of virtual forces and torques, we obtain the theorem of complementary 
virtual work [52j : 


- E • /b 

ijeCs 


■ Fi+uji 


n 


+ E vtrXn. 


(15) 
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FIG. 13. (Color online) Strain scale from collapse of C'(e). 
Lines as in Fig.5a of the main text. 


where /„ and 




^ fij 


are the net virtual contact force and contact torque on particle i, 


and x nij on the boundary. Eq. (15) applies only in between collisions. It is important to 


stress that Eq. ( |15[ ) holds for any virtual force field {fa}- To compute the nonlinearity in flow, we will use the virtual 
work theorem with {dfa/de} as the virtual ‘forces’. This allows us to identically remove the leading order terms in 
flow and consider only those that evolve with strain. The basic equation is then 


E - E ”« %+E 




ijeCs 


dFi - dn ^ 

~r~ ' T ~T~ ' 
de de 


dfta 


E^.ij..E = x/L 


(16) 


where contacts are denoted as a. Here the last term is needed to precisely cancel the dfia/de terms that appear in 
dTi/de when expanded. Under constant stress boundary conditions, we can fix all the boundary forces, so that the 
LHS vanishes. The terms on the RHS can be simplihed using the definition of floppy modes, and the fact that at 
sliding contacts we have |/J| = fJ-pfa ■ After a long computation, in d = 2 we can rewrite Eq. (16) exactly as 


o=E 

aeCs 


jfN fN 

a J cx. f-i \2 


-E 


N 


eSr ' 




E 


dF, 

de 


V, 


dfi 

de 


■ Wj 


(17) 


where Ua is the magnitude of sliding velocity at contact a. In d = 3 there are several additional terms that are not 
expected to be important, for example involving the slight difference between sliding directions and the directions of 
tangential forces. 

When flow is along floppy modes, velocities have a characteristic scale dU; we will assume that the angular and 
linear velocities have the same scale, 6 uj ~ SVfD. Then since y < 1, in terms of scaling we have 


0 : 


xNct^pSV^ - Nct 6V^ + 

de e de 


(18) 


where we used that dFi/de-Vi > 0. Under constant stress BCs, the dp/de term is negligible (more precisely, it vanishes 
up to correlations between Ua and df^/de). Then we find 


dF ZPc,r 

Ik “ 


(19) 


as stated in the main text. Eq.(19) indicates how quickly conhgurations flow out of equilibrium along floppy modes. 


and applies for both viscous and inertial dynamics. The magnitude of unbalanced forces, F, can itself be written in 
terms of geometrical quantities, but this depends on the dynamics. 
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